PHREEQCをゼロから始める #8：Python で溶解度データを解析・可視化する

#7 で出力した Gibbsite 溶解度データを pandas で読み込み、matplotlib / plotly でインタラクティブなダイアグラムを生成する。PHREEQC と Python の連携ワークフローを完成させる。
PHREEQC
Python
地球化学
作者

DeepFlow

公開

2026年5月4日

はじめに：なぜ Python と連携するのか

PHREEQC の USER_GRAPH は手軽にグラフを描けるが、出版品質の図・インタラクティブ可視化・統計処理には Python が圧倒的に強い。 #7 で出力した gibbsite_solubility.txt を Python で読み込み、以下を実現する：

🐼
pandas
SELECTED_OUTPUT の
タブ区切りファイルを
スマートに読み込む
📈
matplotlib
出版品質の
溶解度ダイアグラムを
描画する
🔵
plotly
ホバー対応の
インタラクティブ図を
ブログに埋め込む
ノート前提
  • #7 の PHREEQC コードを実行済みで gibbsite_solubility.txt が手元にあること
  • Python 3.9 以上 + pandas, matplotlib, plotly がインストール済みであること
    （pip install pandas matplotlib plotly kaleido ）kaleidoは画像ファイルの出力
  • pipは Python のパッケージ管理ツールで、ターミナル（コマンドプロンプト）で上記コマンドを実行することで必要なライブラリをインストールできる。
  • 例：pip install phreeqpy で phreeqpy をインストールできる。上記のコマンドは pandas、matplotlib、plotly、kaleido を一括でインストールする。

Step 1：SELECTED_OUTPUT ファイルの構造を理解する

PHREEQC の SELECTED_OUTPUT が出力するファイルはタブ区切りのテキストである。 先頭行がヘッダー、2行目以降がデータになる。

# gibbsite_solubility.txt（抜粋）
sim state   soln    dist_x  pH  Al Al+3    AlOH+2  Al(OH)2+    Al(OH)3 Al(OH)4-    Gibbsite
1   i_soln  1   -1  7.000   0.000e+00   ...
2   react   1   -1  3.000   3.000   7.76e-03    -2.110  -5.270  -7.970  -9.710  -13.5   ...
3   react   2   -1  4.000   7.76e-04    -3.110  -5.270  -6.970  -9.710  -13.5   ...
...
重要ヘッダー名の注意点

-activities で指定した化学種は log₁₀（活量） で出力される。 一方 -totals で指定した Al は mol/kgw（線形値） で出力される。 可視化前に np.log10() で変換が必要な列と不要な列を区別しよう。

場合によってはヘッダーが　「pH Al・・・」に始まる場合がある。

Step 2：pandas でデータを読み込む

# ============================================================
#  phreeqc_read.py
#  PHREEQC SELECTED_OUTPUT ファイルの読み込みと整形
# ============================================================
import pandas as pd
import numpy as np
from pathlib import Path

# ---- ファイル読み込み ----
fp = Path("gibbsite_solubility.txt")

df_raw = pd.read_csv(
    fp,
    sep="\t",        # タブ区切り
    comment="#",     # # 行をスキップ
    skipinitialspace=True,
)

print(f"Shape: {df_raw.shape}")
print(df_raw.head())
print("\nColumns:", df_raw.columns.tolist())
> Shape: (13, 10) or (10, 10)
> Columns: ['sim', 'state', 'soln', 'dist_x', 'pH', 'Al', 'la_Al+3',
         'la_AlOH+2', 'la_Al(OH)2+', 'la_Al(OH)3', 'la_Al(OH)4-', 'Gibbsite'] 
# ---- 整形：初期状態行を除き、pHデータ行だけ残す ----
df = (
    df_raw
    #.query("state == 'react'")          # 平衡計算結果だけ
    .rename(columns={
        "Al": "Al_total_mol",  # total Al（線形値）
        "la_Al+3":        "log_Al3",
        "la_AlOH+2":      "log_AlOH2",
        "la_Al(OH)2+":    "log_AlOH2p",
        "la_Al(OH)3":     "log_AlOH3",
        "la_Al(OH)4-":    "log_AlOH4",
    })
    .reset_index(drop=True)
)

# ---- total Al を log スケールに変換 ----
df["log_Al_total"] = np.log10(df["Al_total_mol"].clip(lower=1e-15))

print(df[["pH", "log_Al_total", "log_Al3", "log_AlOH2p", "log_AlOH3","log_AlOH4"]].to_string(index=False))
pH log_Al_total log_Al3 log_AlOH2p log_AlOH3 log_AlOH4
3.12 0.45 -0.98 -5.03 -8.83 -11.52
4.00 -3.71 -3.89 -6.00 -8.83 -10.56
5.00 -6.44 -6.89 -7.00 -8.83 -9.56
6.00 ★ -7.80 -9.89 -8.00 -8.83 -8.56
7.00 -7.52 -12.89 -9.00 -8.83 -7.56
8.00 -6.55 -15.89 -10.00 -8.83 -6.56
9.00 -5.55 -18.89 -11.00 -8.83 -5.56
10.00 -4.55 -21.89 -12.00 -8.83 -4.56
11.00 -3.54 -24.89 -13.00 -8.83 -3.56
12.00 -2.50 -27.89 -14.00 -8.83 -2.56

Step 3：matplotlib で出版品質の図を描く

Step 2の下に以下のコードをコピー＆ペーストすること。

# ============================================================
#  phreeqc_plot_matplotlib.py
# ============================================================
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import numpy as np
import pandas as pd

# ---- フォント設定（日本語環境） ----
plt.rcParams.update({
    "font.family": "sans-serif",
    "font.sans-serif": ["MS Gothic", "Noto Sans CJK JP", "IPAexGothic", "DejaVu Sans"],
    "axes.unicode_minus": False,
    "figure.dpi":      150,
})

# ---- データ（dfはStep2で作成済み） ----
ph = df["pH"].values

species = {
    "Total Al":      ("log_Al_total", "#111827", 2.8, "-",  "o"),
    r"Al$^{3+}$":    ("log_Al3",      "#DC2626", 1.8, "-",  "s"),
    r"AlOH$^{2+}$":  ("log_AlOH2",   "#EA580C", 1.6, "--", "^"),
    r"Al(OH)$_2^+$": ("log_AlOH2p",  "#CA8A04", 1.6, "-.", "D"),
    r"Al(OH)$_3$(aq)":("log_AlOH3",  "#16A34A", 1.6, ":",  ""),
    r"Al(OH)$_4^-$": ("log_AlOH4",   "#2563EB", 1.8, "-",  "v"),
}

fig, ax = plt.subplots(figsize=(9, 6))

# ---- Gibbsite 安定域の背景色 ----
ax.axvspan(5.5, 8.5, alpha=0.07, color="#16A34A", label="_nolegend_")
ax.text(7.0, -2.3, "Gibbsite\n安定域", ha="center", va="top",
        fontsize=9, color="#15803D", style="italic")

# ---- 各化学種ライン ----
for label, (col, color, lw, ls, marker) in species.items():
    y = df[col].replace([np.inf, -np.inf], np.nan).values
    valid = ~np.isnan(y) & (y > -12)          # 検出限界以上のみ
    if marker:
        ax.plot(ph[valid], y[valid], color=color, lw=lw,
                ls=ls, marker=marker, ms=5, label=label)
    else:
        ax.plot(ph[valid], y[valid], color=color, lw=lw,
                ls=ls, label=label)

# ---- 最小値アノテーション ----
min_idx = df["log_Al_total"].idxmin()
min_ph  = df.loc[min_idx, "pH"]
min_log = df.loc[min_idx, "log_Al_total"]
ax.annotate(
    f"最小値\npH ≈ {min_ph:.1f}\n10⁻⁷ mol/L",
    xy=(min_ph, min_log),
    xytext=(min_ph + 1.5, min_log + 1.2),
    fontsize=9, color="#92400E",
    arrowprops=dict(arrowstyle="->", color="#D97706", lw=1.5),
    bbox=dict(boxstyle="round,pad=0.3", fc="#FEF3C7", ec="#D97706", alpha=0.9),
)

# ---- 軸・装飾 ----
ax.set_xlim(3, 14)
ax.set_ylim(-11, 0.5)
ax.set_xlabel("pH", fontsize=13)
ax.set_ylabel(r"$\log$ [Al]  (mol/kgw)", fontsize=13)
ax.set_title("Gibbsite 溶解度ダイアグラム（Al–H₂O系, 25°C）", fontsize=14, pad=12)
ax.xaxis.set_major_locator(ticker.MultipleLocator(1))
ax.yaxis.set_major_locator(ticker.MultipleLocator(1))
ax.grid(True, which="major", ls="--", lw=0.5, color="#E5E7EB")
ax.legend(loc="upper right", fontsize=10, framealpha=0.95)

plt.tight_layout()
plt.savefig("gibbsite_solubility.png", dpi=300, bbox_inches="tight")
plt.savefig("gibbsite_solubility.svg", bbox_inches="tight")   # ベクター形式
plt.show()
print("✅ 保存完了：gibbsite_solubility.png / .svg")

Figure of the Result of Step 3

Step 4：優勢化学種マップ（speciation diagram）

溶解度だけでなく、どの化学種が優勢かを色で示す図も地球化学では頻出する。

# ============================================================
#  phreeqc_speciation_map.py
#  各 pH での優勢 Al 化学種を色でマッピング
# ============================================================
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import numpy as np
import pandas as pd

# 各 pH で最大活量の化学種を判定
activity_cols = ["log_Al3", "log_AlOH2", "log_AlOH2p", "log_AlOH3", "log_AlOH4"]
species_names  = ["Al³⁺", "AlOH²⁺", "Al(OH)₂⁺", "Al(OH)₃(aq)", "Al(OH)₄⁻"]
species_colors = ["#DC2626", "#EA580C", "#CA8A04", "#16A34A", "#2563EB"]

df_act = df[activity_cols].replace([np.inf, -np.inf], np.nan).fillna(-99)
dominant_idx = df_act.values.argmax(axis=1)

fig, ax = plt.subplots(figsize=(10, 2.2))

for i, row in df.iterrows():
    ph_val = row["pH"]
    d_idx  = dominant_idx[i]
    color  = species_colors[d_idx]
    rect   = mpatches.FancyBboxPatch(
        (ph_val - 0.45, 0.1), 0.9, 0.8,
        boxstyle="round,pad=0.05",
        facecolor=color, edgecolor="white", linewidth=1.5, alpha=0.85,
    )
    ax.add_patch(rect)
    ax.text(ph_val, 0.5, species_names[d_idx],
            ha="center", va="center", fontsize=8, color="white", fontweight="bold")

ax.set_xlim(2.5, 12.5)
ax.set_ylim(0, 1)
ax.set_xlabel("pH", fontsize=12)
ax.set_yticks([])
ax.set_xticks(range(3, 13))
ax.set_title("Al 優勢化学種マップ（Gibbsite 平衡, 25°C）", fontsize=12, pad=8)

# 凡例
patches = [mpatches.Patch(color=c, label=n)
           for c, n in zip(species_colors, species_names)]
ax.legend(handles=patches, loc="upper center",
          bbox_to_anchor=(0.5, -0.35), ncol=5, fontsize=9, framealpha=0.9)

plt.tight_layout()
plt.savefig("al_speciation_map.svg", bbox_inches="tight")
plt.show()

計算結果に基づくAl優勢化学種マップ シミュレーション結果により、pH 6付近での Al(OH)2+の出現や、広範なアルカリ域での Al(OH)4の優勢状況が確認できます。

Step 5：完全な解析ワークフロー（one-file スクリプト）

上記 Step 1〜5 をひとつのスクリプトにまとめた完全版で、このコードをコピー＆ペーストして図を確認しよう。#7からの成果物である。

#  gibbsite_full_analysis.py
#  PHREEQC #7 結果の完全解析スクリプト
#  使い方: python gibbsite_full_analysis.py gibbsite_solubility.txt
# ============================================================
import sys
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import matplotlib.patches as mpatches

# ---- フォント設定（日本語環境） ----
plt.rcParams.update({
    "font.family": "sans-serif",
    "font.sans-serif": ["MS Gothic", "Noto Sans CJK JP", "IPAexGothic", "DejaVu Sans"],
    "axes.unicode_minus": False,
    "figure.dpi":      150,
})

# ---- 定数 ----
SPECIES = {
    "Total Al":       ("log_Al_total", "#111827", 2.8, "-"),
    r"Al$^{3+}$":     ("log_Al3",      "#DC2626", 1.8, "-"),
    r"AlOH$^{2+}$":   ("log_AlOH2",   "#EA580C", 1.6, "--"),
    r"Al(OH)$_2^+$":  ("log_AlOH2p",  "#CA8A04", 1.6, "-."),
    r"Al(OH)$_3$":    ("log_AlOH3",   "#16A34A", 1.6, ":"),
    r"Al(OH)$_4^-$":  ("log_AlOH4",   "#2563EB", 1.8, "-"),
}

ACT_COLS  = ["log_Al3", "log_AlOH2", "log_AlOH2p", "log_AlOH3", "log_AlOH4"]
SP_NAMES  = [r"Al$^{3+}$", r"AlOH$^{2+}$", r"Al(OH)$_2^+$",
             r"Al(OH)$_3$(aq)", r"Al(OH)$_4^-$"]
SP_COLORS = ["#DC2626", "#EA580C", "#CA8A04", "#16A34A", "#2563EB"]

def load_data(filepath: str) -> pd.DataFrame:
    """SELECTED_OUTPUT を読み込み整形する"""
    df = pd.read_csv(filepath, sep="\t", comment="#", skipinitialspace=True)
    #df = df.query("state == 'react'").reset_index(drop=True)
    df = df.rename(columns={
        "Al": "Al_total_mol",
        "la_Al+3": "log_Al3", "la_AlOH+2": "log_AlOH2",
        "la_Al(OH)2+": "log_AlOH2p", "la_Al(OH)3": "log_AlOH3", "la_Al(OH)4-": "log_AlOH4",
    })
    df["log_Al_total"] = np.log10(df["Al_total_mol"].clip(lower=1e-15))
    return df

def plot_solubility(df: pd.DataFrame, save: bool = True):
    """溶解度ダイアグラムを描画"""
    fig, ax = plt.subplots(figsize=(9, 6))
    ax.axvspan(5.5, 8.5, alpha=0.07, color="#16A34A")
    ax.text(7.0, -2.0, "Gibbsite\n安定域", ha="center",
            fontsize=9, color="#15803D", style="italic")

    for label, (col, color, lw, ls) in SPECIES.items():
        y = df[col].replace([np.inf, -np.inf], np.nan).values
        mask = ~np.isnan(y) & (y > -12)
        ax.plot(df["pH"].values[mask], y[mask],
                color=color, lw=lw, ls=ls, label=label, marker="o", ms=4)

    # 最小値アノテーション
    mi = df["log_Al_total"].idxmin()
    ax.annotate(
        f"最小値 pH≈{df.loc[mi,'pH']:.1f}",
        xy=(df.loc[mi, "pH"], df.loc[mi, "log_Al_total"]),
        xytext=(df.loc[mi, "pH"] + 1.5, df.loc[mi, "log_Al_total"] + 1.5),
        fontsize=9, color="#92400E",
        arrowprops=dict(arrowstyle="->", color="#D97706"),
        bbox=dict(boxstyle="round", fc="#FEF3C7", ec="#D97706"),
    )

    ax.set(xlim=(3, 14), ylim=(-11, 0.5),
           xlabel="pH", ylabel=r"$\log$ [Al] (mol/kgw)",
           title="Gibbsite 溶解度ダイアグラム（Al–H₂O系, 25°C）")
    ax.xaxis.set_major_locator(ticker.MultipleLocator(1))
    ax.yaxis.set_major_locator(ticker.MultipleLocator(1))
    ax.grid(True, ls="--", lw=0.5, color="#E5E7EB")
    ax.legend(fontsize=10, loc="upper right")
    plt.tight_layout()
    if save:
        fig.savefig("gibbsite_solubility.svg", bbox_inches="tight")
        print("✅ gibbsite_solubility.svg を保存")
    plt.show()

def plot_speciation_map(df: pd.DataFrame, save: bool = True):
    """優勢化学種マップを描画"""
    df_act = df[ACT_COLS].replace([np.inf, -np.inf], np.nan).fillna(-99)
    dom = df_act.values.argmax(axis=1)

    fig, ax = plt.subplots(figsize=(10, 2.2))
    for i, row in df.iterrows():
        c = SP_COLORS[dom[i]]
        rect = mpatches.FancyBboxPatch(
            (row["pH"] - 0.45, 0.1), 0.9, 0.8,
            boxstyle="round,pad=0.05",
            facecolor=c, edgecolor="white", lw=1.5, alpha=0.85,
        )
        ax.add_patch(rect)
        ax.text(row["pH"], 0.5, SP_NAMES[dom[i]],
                ha="center", va="center", fontsize=8,
                color="white", fontweight="bold")

    patches = [mpatches.Patch(color=c, label=n)
               for c, n in zip(SP_COLORS, SP_NAMES)]
    ax.set(xlim=(2.5, 12.5), ylim=(0, 1),
           xlabel="pH", title="Al 優勢化学種マップ（Gibbsite 平衡, 25°C）")
    ax.set_yticks([])
    ax.set_xticks(range(3, 13))
    ax.legend(handles=patches, loc="upper center",
              bbox_to_anchor=(0.5, -0.35), ncol=5, fontsize=9)
    plt.tight_layout()
    if save:
        fig.savefig("al_speciation_map.svg", bbox_inches="tight")
        print("✅ al_speciation_map.svg を保存")
    plt.show()

# ---- メイン ----
if __name__ == "__main__":
    fp = sys.argv[1] if len(sys.argv) > 1 else "gibbsite_solubility.txt"
    df = load_data(fp)
    print(df[["pH", "log_Al_total", "log_Al3", "log_AlOH4"]].to_string(index=False))
    plot_solubility(df)
    plot_speciation_map(df)

実行は一行：

python gibbsite_full_analysis.py gibbsite_solubility.txt

まとめ：PHREEQC × Python ワークフロー全体像

PHREEQC .pqi 実行 SELECTED_OUTPUT gibbsite_solubility.txt pandas 読込・整形 matplotlib 静的図 (SVG)
Step 1–2
PHREEQC 実行
→ txt 出力
→ pandas 読込
Step 3
matplotlib で静的図
Step 4–5
優勢種マップ
one-file スクリプトで
自動化完了

次回 #9 では、イオン強度と活量係数 — 海水と真水で計算結果が変わる理由について調べる。

参考文献（References）

Appelo, CAJ, と Dieke Postma. 2005年. Geochemistry, groundwater and pollution. Second. Balkema, Rotterdam, p. 634.
Parkhurst, David L, と CAJ Appelo. 2013年. Description of input and examples for PHREEQC version 3—A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations. US Geological Survey Techniques; Methods, book 6, chap. A43, 497 p.
Yamamoto, S. 1983年. Method of the groundwater survey. Kokon Shoin, Tokyo (in Japanese), 490 p.
Yang, Heejun, T Mishima, S Katazakai, と M Kagabu. 2023年. 「Analytical approach using a chemical equilibrium formula and geochemical modeling for alkalinity measurements of small natural water samples」. Applied Geochemistry 148: 105535.

DeepFlow | 地球科学シミュレーションの深みへ